Leading-Order Auxiliary Field Theory of the Bose-Hubbard Model 
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We discuss the phase diagram of the Bose-Hubbard (BH) model in the leading-order auxiliary field 
(LOAF) theory. LOAF is a conserving non-perturbative approximation that treats on equal footing 
the normal and anomalous density condensates. The mean-field solutions in LOAF correspond to 
first-order and second-order phase transition solutions with two critical temperatures corresponding 
to a vanishing Bose-Einstein condensate, Tc, and a vanishing diatom condensate, T*. The second- 
order phase transition solution predicts the correct order of the transition in continuum Bose gases. 
For either solution, the superfluid state is tied to the presence of the diatom condensate related to 
the anomalous density in the system. In ultracold Bose atomic gases confined on a three-dimensional 
lattice, the critical temperature Tc exhibits a quantum phase transition, where T^ goes to zero at a 
finite coupling. The BH phase diagram in LOAF features a line of first-order transitions ending in 
a critical point beyond which the transition is second order while approaching the quantum phase 
transition. We identify a region where a diatom condensate is expected for temperatures higher 
than Tc and less than To, the critical temperature of the non-interacting system. The LOAF phase 
diagram for the BH model compares qualitatively well with existing experimental data and results 
of ab initio Monte Carlo simulations. 

PACS numbers: 03.75.Hh, 05.30.Jp, 67.85. Be 



I. INTRODUCTION 

The Bose-Hubbard (BH) model has been the subject 
of broad theoretical [1-7] and experimental [8-10] inter- 
est. In addition to being a challenge for many-body theo- 
ries [11[, its realization in ultracold atoms [9[ opened op- 
portunities for studying the BH model in a controllable 
and accurate way. Developing and improving mean-field 
descriptions for the BH model has been an important task 
since Fischer et al. [1[ discussed the zero-temperature 
mean-field phase diagram for the BH model. Studies 
of the BH model have been summarized in many text- 
books [12-14[. One may also test various many-body 
theoretical techniques using the BH model and bench- 
mark those methods. Here we follow this tradition and 
study the BH model using a well-developed theoretical 
framework. 

Recently, we introduced a leading-order auxiliary field 
(LOAF) theory for a homogeneous system of ultracold 
gas of bosonic atoms [15, 16[. To derive this formal- 
ism, we used the Hubbard-Stratonovitch transforma- 
tion [17, 18[ to introduce auxiliary fields related to the 
normal and anomalous density condensates. Path inte- 
gral methods were used to obtain a leading-order expan- 
sion of the partition function using the auxiliary fields 
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to organize the expansion method. The resulting non- 
perturbative mean-field theory produces a conserving 
and gapless approximation that is applicable to large in- 
terval of coupling-constant values, satisfies the Goldstone 
theorem, yields a Bose-Einstein transition that is second 
order, and predicts a positive shift in the critical tem- 
perature, Tc, consistent with other similar methods [19[. 
The relation of the LOAF theory to the Goldstone theo- 
rem and the Higgs mechanism was discussed in Ref. 20. 
The behavior of the LOAF theory near the critical point 
is discussed in Ref. 21 and the relation to superfluidity 
and the Josephson relation is discussed in Ref. 22, where 
we showed that the superfluid density in LOAF is propor- 
tional with the square of the anomalous-density diatom 
condensate. The latter is analogous with the Cooper-pair 
condensate in the BCS mean-fleld theory of dilute Fermi 
gases [23[. 

In this paper we develop the LOAF theory of the Bose- 
Hubbard model, which has been used to study the physics 
of ultracold Bose atoms in optical lattices [4-7[. Per- 
haps the most salient feature of the BH model is the 
prediction of a superfluid to Mott insulator phase tran- 
sition. The latter was demonstrated experimentally in 
three-dimensional optical lattices by Trotzky et al. [10[ 
by observing the suppression of the critical temperature 
for superfluidity near the Mott transition. These ex- 
perimental results were showed to compare nicely with 
theoretical predictions based on quantum Monte Carlo 
simulations [5[. Monte Carlo also predicted the critical 
interaction strength for which the critical temperature 
goes to zero and the phase transition is purely quantum 
in character [5[. 



We will show that below the critical temperature, 
Tc, where the usual Bose-Einstein condensate vanishes. 
LOAF is identical with the Hamiltonian version of LOAF 
introduced recently by Kleinert, Narzikulov, and Rakhi- 
mov [24] and referred to as the "two-collective" field the- 
ory for the BH model. Above Tc, LOAF has two pos- 
sible solutions corresponding to first-order and second- 
order phase transitions, respectively. For both solu- 
tions, the superfluid state is indicated by the presence 
of an anomalous-density diatom condensate in the sys- 
tem. In ultracold Bose atomic gases confined on a three- 
dimensional lattice, LOAF predicts that the critical tem- 
perature Tc exhibits a quantum phase transition (QPT), 
where Tc goes to zero at a finite coupling. The BH phase 
diagram in LOAF features a line of first-order transitions 
ending in a critical point beyond which the transition 
is second order while approaching the QPT limit. We 
identify a region where a diatom condensate is expected 
for temperatures Tc < T < Tq. Here, Tg is the critical 
temperature of the non-interacting system. Overall, the 
LOAF phase diagram is very similar to the superfluid 
to Mott insulator transition in systems of ultracold Bose 
atoms trapped in optical lattices and compares qualita- 
tively well with results of ab initio Monte Carlo simula- 
tions [5] and available experimental data [10]. We ana- 
lyze numerically the properties of the LOAF theory in 
the weak-coupling limit. 



II. THE BOSE-HUBBARD MODEL 



energy between atoms. The tight-binding approximation 
assumes that the field (/>(x, t) can be expanded in nor- 
malized eigenfunctions ipipc) of an atom trapped in the 
lattice at positions Xi — a i. 



(x,t)=^0i(t)^(x~Xi) 



(3) 



where ip{x) satisfies 
2m 



+ V{x) V(x)=^V'(x) 



(4) 



where E is the ground state energy of a trapped atom of 
mass m. Here i = {ix,iyiiz) s^re triplets of integers, each 
running from 1 to A^^. Inversion of (3) gives 



(j)i{t)^ / d^a;V*(x-Xi)0(x,i). 



(5) 



So using expansion (3) , and keeping overlaps with nearest 
neighbors, we find 



A^x<j)*{^,t) 



2m 

|2 



y(x) 



(x,i) 



(6) 



i 



We consider the case of N bosonic atoms trapped in a 
three-dimensional cubic lattice. 



where k = (1,0, 0), (0, 1,0), (0, 0, 1) is the displacement 
by one unit in the {x, y, 2:)-directions, and J is the overlap 
integral. 



Real time formulation 



The path integral for the boson field (j>{x, t) is given by 



z[].r 



iW 



[].3']/h ^ ffB(f>B<t)*e'^^'f'''l'''^-^^'^/^, 



S[' 



;.?,j 



dtLl 



\J,J 



(1) 



where the Lagrangian is 
ih 



(2) 



d^x { 0* (x, t) [dtH^, t)] - [dtcb* (x, t)] 0(x, t) } 



d3x0*(x,t) 



2y72 



2m 



l/(x) 0(x,i) 



J 



d^xV'(x + a)y(x)V'(x) > 0, 



(7) 



The particle-particle interaction is assumed to be a short 
range contact interaction. 



1 



d^x d^x' (j)* (x, t) (P* (x', t) t/(x, x') 0(x', t) (/)(x, t) 



U 



Ei-^iW 



(8) 



So in the tight-binding approximation, the Lagrangian is 
given by 



ih 



(9) 



d'x d^x' r (x, t) r (x', t) C/(x, x') </>(x', t) 0(x, t) = y E { '^i (^) t^* '^i (^)] " t^* "^i (^)] ^^ (*) } 



+ y d3x[j*(x,i)0(x,t)-|-j(x,i)r(x,t)] . 

Variation of the action yields a Schrodinger equation for 
the field 0(x, t). Here V{x.) is the parodic potential cre- 
ated by the optical lattice and C/(x, x') is the interaction 



-E{s i0i(t)i' - ^EK*(^) <^i+-w + '^i+KW -^i w] } 



u 



Ei'^'Wi'+Et-^iW'^iW+jiW'^iW]' 



which we can write as: 



III. LOAF FORMALISM 



ih 



(10) 



= TE{'^iW[^*'^iW]-[^*'^iW]<^iW} 



u 



i i 

i.K. 
i 

Here d is the number of spatial dimensions. The constant 
energy term proportional to i? — 2dJ can be eliminated 
by changing variables to 



,/,i(t)=e-«(^-2'i^)*/''0i(t), 



(11) 



which simply changes the energy scale. Then in terms of 
these new variables, (10) becomes 






(12) 



u 



Ei^'Wi'+E[j'iW^iW+j'iW<^iW: 



Eq. (12) is the usual form of the Bose-Hubbard La- 
grangian. From now on we drop the tilde notation. 



B. Imaginary time formulation 

The imaginary time action is obtained from the real 
time action by the mapping t i— > —ihr and L i— )• — Le. 
The partition function Z for the BH model is then written 
as 



*i - ^-/3f^[i,r] _ 



Z[j,.f] = e 



T)(l)B(t>*e 



-SE[<f',r;j,r] 



SE[c^,r-:J,J*]^ / dTLE[0,0*;j,j*], (13) 
Jo 

where the Euclidean Lagrangian is 

iE[0>*;j,r] (14) 

= E { ^ { -^i (^) ^^rMr)] - [dr^l^Ur)] Mr) } 

i 

+ J^{2|0i(Ol'-[0i(r)0i+«(r) + 0;^+,(r)0i(r)]} 

K 

+ I \Mr) I' -Ml0i(r)|^ - jT(r)0i(r) - Ji iT)<p;ir) } . 

Here we have dropped the tilde notation and introduced 
a chemical potential fj,. 



In the leading-order auxiliary field (LOAF) method, we 
introduce two auxiliary fields Xi(''') ^iid Ai(T) by means 
of the Hubbard-Stratonovitch transformation [17, 18]. In 
our case, the auxiliary-field Lagrangian density takes the 
form 

Lauxi*, A] - E { ^ I ^-M -U'l^fir) r (15) 

i 

-l^[xiir)-UV2\Mr)\'T}, 

which we add to Eq. (14). We show in Ref. 16 that this 
choice, in the weak coupling limit, agrees with Bogoli- 
ubov theory [25, 26]. The action is then becomes 



= \ ly j^'Yl ^Ur)gzlir.T')'^^{T') 



(16) 



dr 



n 



^''^^''-"% ^(^)^,M 



+ K; 



jY.[2\Ut)?~[m)k+M+'^UMk{t)]} with 



"/(^)Ai(r)} 



QVliry) 



(17) 



5{t,t') 






where 



/iij-JVij+5ij[N/2x(r)-M], 
ViJ=E{2'5iJ-['5iJ+« + '5i+«j]} 



(18a) 
(18b) 



Here we have introduced currents which we write as Ji(T) 
and Ki{T) and a notation. 



$i(r) 



Ur) 



Mr) 



Mr) 



,<J>Ur)J ' "''-'-KjUr), 

for the particle fields and currents, and a notation 



(19) 



A,{r)\ h{T) 

Ai(T) = Xi{r) , Mr) = fcoi(r) | , (20) 

VUr)) XKir), 

for the auxiliary fields and currents. The generating func- 
tional for the fields is written as a path integral over all 
the fields 

Z[JJ<] = e-Z'f^I^.^l = /'/'D$DAe-^-[*'^'-^'^l . (21) 

The action is now quadratic in the 0i(T) which can be 
integrated out, giving 



Z[J,K]= /"DAe-^-^"!^^'^'^! 



(22) 



where now 



5cff[A;J,X] 

\ lfdTdr'Y,jl{T)g,:,{Ty)j^{T') 



(23) 



2jJ0 



l>U'^ 



xi{r)-\A,{r)Y 



2U 



+ KI{t)A,{t) 



lTr[ln[gr.i(r,r)]]} 



Expanding the effective action, 



dr 



n 



6S,s[^;J,K] 



SAiir) 



(24) 
(Ai(r)-Ai(r)) 



^/i--S[g 



S^S,s[A;J,K] 



2 77o"'"' ^ L 5Ai(a;)(5Aj(a;') 

x(Ai(r)-Ai(r))(Aj(T')-Aj(r')) + 
about the stationary points A = A, defined by 

SAiir) J A 

and computing the remaining path integral by the 
method of steepest descent, we find 



(25) 



r/3 



(26) 



I liy d^' E ^i (^) ^i j' (^' ^') ^J (^' 

Vr X?(r)-|^i(r )p , ,,t 
^ I 2C/ 



/3 

dr 



+ i^J(T)Ai(r) 



iTr[hr[(?r,i(r,r)]]-^Tr[ln[2?rji(r,r)]]} 



where 



^M(r,r') = 



d^S,s[A;J,K] 



5Ai(T)5Aj(r') 
Here $i(T) is defined as the solution of 



hir)^ /dr'^t;,j(r,r')Jj(r'), 

Jo : 



(27) 



(28) 



and is a functional of the currents {J,K). Explicitly, 
the stationary points are defined by the solutions of the 
equations 

^ = ^|.t(r)$,(r) + iTr[gu(r,r)]+fcoi(r), 



Aiir) 1 



2U 



2 $I(T)a+$i(r) + -Tr[f7_gu(r,r)] - fci(T) , 



(29a) 

(r), 
(29b) 



and are functionals of the currents and a± are the Pauli 
matrices. Introducing the Legendre transformation, 

/3V;ff[$,A]= /drV[j/(r)$i(r) + if/(r)Ai(T)] 
Jo ■ 

-f5n[J,K], (30) 

we obtain the thermodynamic effective potential, 

VM^,A] = ^/^drdr'X; $I(r)e.-;ji(r,r') $j(r') 



+ 



U^U 



2U 
lTr[ln[t;-.i(r,r)]]}, 



(31) 



where we have dropped the trace-log term involved the V 
propagator, which is higher order in our expansion. The 
currents are now given by derivatives of yoff['I', A] with 
respect to the fields. 



^'(^) = ^^^^i;Fr' ^^'(^) = ^^ 



914ff[^,A] 
dAiir) 



The thermodynamic potential is evaluated at zero cur- 
rents, which is at the minimum of V^fi[<i>, A]. The aver- 
age particle number is given by 



N = - 



dVos[^,A] 



(32) 



evaluated at the minimum of the effective potential. 



IV. HOMOGENEOUS SYSTEMS 

For homogeneous lattice systems in equilibrium, the 
fields are independent of r and i. Expanding the inverse 
Green function in a three dimensional Fourier series, 

C-^(t t'] - ^ V Q-^ j[27rk.(i-j)/W,-c^„(r-r')] 
^ k,n 

(33) 
where w„ — 2TTn//3 are the Bose Matsubara frequencies. 
Here k = (fc^, ky, kz) is a triplet of integers, each running 
from —Ns/2 to Ns/2 — 1. The total number of sites in 
the cubic box is N^ and the filling factor, i/, is defined 
to be the number of particles per site, f = N/N^. From 
Eq. (6), the Fourier transform of the Green function is 
given by 



—' Ic 71 



k.n 



fk + X -^^n 



-A 



-A* Ck + x' + "i^nj 



(34) 



where we have put x' = v 2 x ~ Mi ^^id the kinetic energy 
is written in terms of the lattice momentum, k, as 

ek = Jk2 = 2J Y^ [l-cos{27rk,/Ns)], (35) 

s—x,y,z 



Using standard techniques, 

r 

2/3 






+ ^t 



Si 



where 



Wfc 



Wfe 



/? 



ln[l 



-/9wfc 1 I 



= v/(£fe + x')'-|^l' 



(36) 



(37) 



The effective potential (3f) for the homogeneous case 
then becomes 



Feff[a>,A]/iV3 



A 



*^2 - 



(38) 

(xi+M)! , HJ! 

4t/ 2C/ 



iaEl^t-^-^^-^'l + ^^tl-^"'"^]}- 



Here we have renormaHzed the effective potential by sub- 
tracting the zero-point energy. The coupling constant is 
finite and does not need to be renormalized. Minimizing 
the effective potential with respect to the fields gives 



{x'~A*)cf, = 0. 



X' + M 
2U 

A 

U ^ 



= \<i>? 



A 

m 



m^i 



2ujk 



E 



[2nk + l] 
2ujk 



(39a) 

(39b) 
(39c) 



where nj. = l/[e^'^'' — f ] , and with the filling factor given 

by 



N 



1 dVM'^,A] 



d^ 



X' + M 
2U 



(40) 



Because of the U(l) invariance of the Lagrangian, at the 
minimum of the potential we can choose to be real. 
Then, from (39a) A is also real since x' is real. We inter- 
pret |(/)p as the number of condensed particles per site, 
and put 



101' 



1^0 






vno 



(41) 



with the condensate fraction, uq = Nq/N. The sums 
over k then omit the k = mode. The gap equations 
then become 



u ~ vna + 



A A 

[7 = ""° + ]v3 



^E'{^[-^ + i]4}^ (-^) 



E 



, [2nfc-H; 
2a; J, 



(42b) 



where Wfc is given in Eq. (37) with e^ given by Eq. (35). 
For comparison purposes, we note that the LOAF ef- 
fective potential per unit volume for the continuum case 

is: 



Kff[1>,A]/y 

|2 



(43) 



,„r-l[A*- + .4Vl-<^ + ^ 



d^k r 1 r 



(2 



Wfc - £*; - X + 



2\ 



ilnll 



-H'^kA _ 



which gives the equations 



(x'-^*)' 

X' + M I 



0, p 

d 

J2n) 



X' + M 
2A ' 
d^k r Sk+x' 
7r)3 I 2wfc 



(44a) 



[2nfe + l]- 



l4} 



A 
A 



f d^k ( 1 , , 1 ^ 

y(2^{2^[2-'= + ^]-2^} 



(44b) 
(44c) 



where the kinetic energy is the usual ek — h^k^ /{2m). 
The differences between the two theories reduce to the 
naive substitution of the integral with the sum over the 
allowed momenta, an extra term in the renormalization 
of the effective potential, and the kinetic energy modifi- 
cation on the lattice. The continuum coupling constant, 
A — ATrfi^ ao/m, corresponds to the Hubbard parameter, 
U, on the lattice. Here uq is the s-wave scattering length 
in the dilute atomic Bose gas. The solutions II(i) and 
Il(ii) correspond to first-order and second-order phase 
transitions, respectively. 



V. RESULTS AND DISCUSSIONS 

The numerical analysis of the solutions space for 
Eqs. (42), leads to three distinct regions in the Bose- 
Hubbard model phase diagram: 

I. The broken symmetry case where 7^ and x' — 
A. Then cj = ^efc(efc + 2%'). In this region, we 
solve the equations [27]: 



v = vnQ + 



ATsE I 



fk+x' 
2ujk 



2nk 



l]-2 



u = """ + m 



X' v:^,[2nfe + i: 



E 



2ujk 



2 } ' (^^^) 
(45b) 



II. The case when = so that uq = 0, and either 
(i) A = so that uj^ = Cfe + x' £^iid 

--4E'-^- (46) 



This solution corresponds to a first-order 
phase transition. Eq. (46) does not depend on 
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FIG. 1. (Color online) Temperature dependence of the con- 
densate fraction, ung, scaled auxiliary fields, x' a.nd A, scaled 
chemical potential, /i, and effective potential, Veft, for a Hub- 
bard interaction parameter value, U/J = f 0, at unity filling, 
V = 1. Here, the temperature, T, is scaled by To, the criti- 
cal temperature of the ideal (non-interacting) Bose-Hubbard 
model. Solutions II(i) and H(ii), corresponds to first- and 
second-oreder phase transitions, respectively. 



the interaction strength and applies for tem- 
peratures, T > Tc, where Tc is the critical 
temperature defined by the zero condensate 
faction limit, uq — > 0, in Eqs. (45). 



(ii) or < v4 < x' so that iOk = ^(^^rTxO^^^^, 
and 



m ^ i 



2ujk 



1 

U 



* k 



2ujk 



2ni 



1] 



1 



1-^} 



(47a) 
(47b) 



This solution corresponds to a second-order 
phase transition. 

III. the normal case where = and A — Q. In this 
case we solve Eq. (46) as in case II(i) above. 



FIG. 2. (Color online) LOAF scaled critical temperatures, Tc 
and T* , and the corresponding scaled normal-density aux- 
ihary fields, x'c = x'dn)^ x'*, and x'c(i), and effective po- 
tentials, as a function of the Bose-Hubbard model coupling 
constant, {U/J), at unity filling, u = 1. We note that only 
the critical temperature, Tc, features a downturn with the 
interaction strength and leads to a quantum phase transi- 
tion (QPT) for {U/J)c ~ 56.07 where Tc goes to zero. In 
addition, the normal-density auxiliary field x!c(i) fo'' solu- 
tion II(i) is not defined for Tc < To, where To ~ 5.59 is 
the critical temperature of the non-interacting Bose system. 
Therefore, LOAF predicts a critical point (CP) at coordi- 
nates Tcp ~ To and (t//J)cp = 46.02. For coupling values 
{U/J)cp < {U/J) < {U/J)c, the transition is second-order 
and LOAF predicts a diatom condensate A 7^ in the absence 
of the usual Bose-Einstein condensate fraction for tempera- 
tures Tc < T < Tq. Because T* > To, the system is in the 
normal phase for T > To in this region. For coupling con- 
stants {U/J) < {U/J)cp the transition is first-order because 
Kff,c(i) < Kff,c(ii) in that region. 



We note that the LOAF solutions for the cubic lattice 
are identical with the LOAF solutions for the continuum 
system [15, 16, 20]. 

To make contact with Ref 24, we convert the finite 
sums over k to integrals by defining q — 2'k/Ns, so that 




10 



FIG. 3. (Color online) Interaction strength dependence of 
the zero-temperature condensate fraction, vtiq, and the cor- 
responding scaled normal-density auxiliary field, XO) a-t unity 
filling, V = 1. 



formally in Eqs. (42) we substitute 



m 



1 ^ 



d^ 



(48) 



This substitution is exact in the limit Ng — )■ oo. In the 
following, all quantities other than the filling factor, v, 
and the condensate fraction, no, can be scaled by J with- 
out loss of generality. 

We define the critical temperature Tc as the point in re- 
gion I where the usual condensate fraction, riQ, vanishes. 
The second critical temperature, T* , is the temperature 
where the diatom condensate. A, vanishes in region Il(ii). 
In the non-interacting limit, [/ — 7> 0, we have T* — > Tc- 

The critical temperature Tc and fields x^ — Ac are 
given by the solution of Eqs. (46) in region I when no = 0: 



^ I 2w„ 



[2n, + 1 



1 

U 



3 y^itq 

a q 



[2n, + l] 



2w„ 



]-l}, (49a) 

, (49b) 



^,/Tc _ 1 



with ujg = y/eq{eq + 2x'c) [28]. From Eqs. (49), the criti- 
cal value of the Hubbard parameter, {U/J)c, is obtained 
by taking the limit Tc — ?► 0. For ly — 1 we obtain the 
critical Hubbard parameter value {U/J)c ~ 56.076, to 
be compared with the critical value of 29.34(2) obtained 
by ab initio quantum Monte Carlo simulations [5]. 
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FIG. 4. (Color online) Comparison of the coupling constant 
dependence of the LOAF critical temperature, Tc, at unity fill- 
ing, V = 1, with experimental [10] and quantum Monte Carlo 
(QMC) results ]5]. The LOAF value of the critical Hubbard 
parameter value, {U/J)c = 56.076, should be compared to the 
QMC critical value, {U/J)c = 29.34(2), reported in Ref. 5. 
LOAF also predicts a critical point at {U/J)cp ~ 46.02. 
The solid and dashed lines indicate first- and second-order 
phase transitions predicted by LOAF theory, respectively. 
The shaded area depicts the region where a diatom conden- 
sate without the usual Bose-Einstein condensate is expected. 
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FIG. 5. (Color online) Weak-coupling limit behavior of the 
critical temperature, Tc, and corresponding critical value of 
the scaled normal-density auxiliary field, Xc, as a function 
of the Hubbard interaction parameter value, U/J, at unity 
filling, v = 1. 
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FIG. 6. (Color online) Ideal gas critical temperature, To, 
critical value of the scaled Hubbard parameter, {U / J)c, and 
the corresponding critical value of the scaled normal-density 
auxiliary field, Xc, as a function of the filling factor, v. 



The critical temperature T* and field %'* are defined 
by the solution of Eqs. (47) when no = and A ~ {): 



V = I I I A^q riq , Uq ^ 



1 



•^JT* „ I ' 



1 

u 



2n„ + 1 " 



d^gi^^ 
2w, 



(50a) 
(50b) 



with UJq = fg+X'*- 

In Fig. 1 we illustrate the temperature dependence of 
the condensate fraction, i^hq, scaled auxiliary fields, x' 
and A, scaled chemical potential, /i, and effective po- 
tential, T4ff, for a Hubbard interaction parameter value, 
U/J = 10, at unity filling, i/ = 1. We find that in 
region II, we have two possible solutions of the LOAF 
equations, as discussed above. The solution II(i) gives 
rise to discontinuities in the temperature dependence of 
the auxiliary fields and chemical potential that lead to 
a discontinuity in the effective potential as well. This 
behavior is characteristic to a first-order phase transi- 
tion. In contrast, the solution Il(ii) corresponds to a 
second-order phase transition, because the temperature 
dependence of %', A, fj,, and VcS is smooth across Tc- 

For convenience, we denote by Xcd) ^^^ value of the 
normal-density auxiliary field corresponding to the first- 
order phase transition solution II(i) (see Eq. (46)) for 
T = Tc, and we introduce the notation Xcia) ~ x'c to 
indicate the value of x' at Tc corresponding to the second- 
order phase transition solution Il(ii). We have, x'cd) ~^ 



x'c(ii) ill the non-interacting limit, [/ ^- 0. We recall 
that Eq. (46) is independent of the Hubbard parameter 
U and is restricted to temperatures T > Tc- In the non- 
interacting limit we have T^, — ^ Tg, where Tq is the critical 
temperature of the non-interacting lattice Bose system, 
and the non-interacting limit corresponds to x' — ^ 0. So, 
we find that Xc(i) -J> in the Hmit {U -> 0, T^ -> Tq}. 
Therefore a first-order phase transition may occur at Tc 
only for an interaction strength, U, that gives a critical 
temperature Tc > Tq. Furthermore, the solution II(i) is 
only possible for a temperature T > Tq. We will use this 
important observation next. 

The coupling constant dependence of the LOAF crit- 
ical temperatures, Tc and T*, the corresponding scaled 
normal-density auxihary fields, Xc(i), x'c{n)=x'c, and x'*, 
and the corresponding effective potentials at Tc are de- 
picted in Fig. 2. We find that the critical temperature T* 
increases monotonically with the BH model coupling con- 
stant, (U/J), whereas the critical temperature, Tc, in- 
creases with the interaction strength for U/J < 10 and 
then decreases with [U/J). At unity filling, v = \, the 
critical temperature, Tc, goes to zero, for a critical value, 
[U / J)c ~ 56.076. It is important to note that in the con- 
tinuum case of a homogenous system of ultracold Bose 
atomic gases neither of the two critical temperatures Tc 
and T* goes to zero in LOAF [21]. It appears that in 
LOAF the presence of a quantum phase transition is re- 
lated to the reduction in the allowed momentum-vector 
phase space. The latter is a consequence of the bosonic 
atoms being spatially confined on the lattice. 

For completeness in Fig. 3 we illustrate the scaled 
Hubbard parameter, U/J, dependence of the zero- 
temperature condensate fraction, vhq, and the corre- 
sponding scaled normal-density auxiliary field, Xo; ^^ 
unity filling, 1^=1. 

As discussed above, the normal-density auxiliary field 
x'c(i) foi' the first-order phase transition solution II(i) is 
not defined for Tc < Tq. That leads to the possibility 
of a critical point (CP) at coordinates Tcp = Tq and 
{U/J)cp ~ 46.02. 

Recalling that in LOAF the superfiuid density is pro- 
portional to the square of the anomalous-density auxil- 
iary field, A (see discussion in Ref. 22), it follows that 
the phase diagram of the Bose-Hubbard model in LOAF 
features two regions: First, for coupling values (U/J) < 
{U/J)cp, both solutions II(i) and Il(ii) are possible and 
we may have either a first-order or a second-order phase 
transition solution. Because T4ff.c(i) < K:ff,c(ii), LOAF 
predicts a first-order phase transition from the super- 
fluid to the normal phase in this region. This scenario 
corresponds to the solution II (i). Second, for coupling 
values {U/J)cp < (U/J) < {U/J)c, we have Tc < Tq 
and the solution II(i) is not possible for temperatures 
Tc < T < Tq. Hence, the transition is second-order as 
described by solution Il(ii). Because T* > Tq for all 
couplings, solution Il(ii) applies only for temperatures 
Tc < T < Tq. In this temperature range LOAF predicts 
a diatom condensate A ^ in the absence of the usual 
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Bose-Einstein condensate fraction and the system is in 
a superfluid state for all temperatures < T < Tq. As 
seen in Fig. 1, for all temperatures Tc < T < T* we have 
Kff,(i) < ^eff,(ii)- Therefore at To the system undergoes 
a first-order phase transition from the superfluid to the 
normal phase. 

The LOAF phase diagram for the BH model can 
be compared with predictions of quantum Monte Carlo 
(QMC) simulations [5]. The depression in the critical 
temperature was observed experimentally in ultracold 
Bose atom systems in three-dimensional lattices [10] and 
the QMC compares well with experiments for couplings 
([// J) < 20. The QPT predicted by Monte Carlo occurs 
for {U/J)c = 29.34(2), so about half the critical value 
predicted by LOAF. In Fig. 4 we show that our LOAF 
results compare qualitatively well with existing experi- 
mental and ab initio QMC results [5, 10]. The shaded 
area in Fig. 4 is the region where a diatom condensate 
is expected to be present in the system in the absence of 
the usual Bose-Einstein condensate ]29]. 

Just like in the continuum case, the LOAF results for 
the lattice show that the critical temperature, T^ departs 
from the ideal gas result (see Fig. 1). Numerical results 
depicted in Fig. 5 show that in the weak-coupling limit 
(ATc) increases linearly with the coupling (U/J) with 
a slope parameter ~ 0.0737. Furthermore, we find that 
the critical value of the auxiliary field, Xc is proportional 
with the square of the coupling, [U / J)"^ . 

Finally, in Fig. 6 we plot the ideal gas critical temper- 
ature. To J critical value of the scaled Hubbard parame- 
ter, [U I J)c, and the corresponding critical value of the 
scaled normal-density auxiliary field, xj,, as a function of 
the filling factor, v. 



VI. CONCLUSIONS 

In this paper we developed the leading-order auxiliary 
field approximation (LOAF) for the Bose-Hubbard model 
corresponding to a system of Bose atoms confined in a 
three-dimensional lattice. The auxiliary-field formalism 
treats on equal footing condensates associated with the 
normal and anomalous densities. 

For temperatures T < Tc we showed that LOAF is the 
same as the "two-collective" field theory introduced by 
Kleinert et al. ]24]. Here T^ is the temperature where 
the usual Bose-Einstein condensate vanishes. For tem- 
peratures T > Tc, LOAF has two possible solutions cor- 
responding to either a first-order or a second-order phase 
transition. For both solutions, the superfiuid state is in- 
dicated by the presence of an anomalous-density diatom 
condensate in the system. 

The BH phase diagram in LOAF features a line of first- 
order transitions ending in a critical point at Tcp = Tq 
and a finite coupling {U/J)cp- In the first-order phase 
transition solution, the diatom condensate auxiliary field 
A vanishes at T^, so the system evolves from a superfiuid 
to a normal phase. First-order phase transition solutions 



are limited to temperatures T > Tq, where Tq is the 
critical temperature of the non-interacting system. 

Beyond the critical point, the transition is second or- 
der. In the case of the second-order phase transition, 
the diatom condensate auxiliary field. A, goes to zero 
smoothly and vanishes at the critical temperature, T*. 
For Tc < T < T* , the system is still in the superfiuid 
phase, because in LOAF the superfiuid density is pro- 
portional to the square of A, and not to the usual con- 
densate fraction, n^. For T > T* , we have ^4 = and 
the system is in the normal phase. This scenario pro- 
vides for the second-order phase transition known to take 
place in liquid helium and dilute gases of Bose atoms. 
The critical temperature T* does not vanish either in 
the continuum or the lattice cases. In the case of the 
BH model, T* is always greater than Tq, so the system 
never reaches T* , but rather exhibits a first-order phase 
transition at Tq from the superfluid to the normal phase. 
Contrary to the conclusions of Kleinert et al. ]24], for 
coupHngs {U/J)cp < {U/J) < {U/J)c and temperatures 
Tc < T < Tq we have a region where a diatom, conden- 
sate is expected in the absence of the usual Bose-Einstein 
condensate. 

For Bose systems on a lattice the critical temperature 
Tc goes to zero for a flnite value of the Hubbard inter- 
action parameter, {U / J)c, indicating a quantum phase 
transition (QPT) similar to the superfluid to Mott in- 
sulator transition. For continuum systems, Tc does not 
vanish ]21] and LOAF predicts no QPT in the case of 
inflnite Bose matter. So in LOAF the QPT is due to 
the spatial confinement of the Bose atoms on the lattice. 
The LOAF phase diagram of the BH model compares 
qualitatively well with existing experimental and ah ini- 
tio quantum Monte Carolo results ]5, 10]. 

It is clear that in order to understand the LOAF phase 
diagram of strongly interacting systems of particles it is 
necessary to extend this non-perturbative theory beyond 
the mean-field level of approximation discussed so far. 
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Appendix A: Non-interacting case 

The non-interacting case, U/J — 0, is recovered by 
setting X — > and A — > in Eq. (38). Then x' — —Mi 
'^fe — ^ Cfe — M- E^iid the effective potential becomes 



Feff[z,T] = i^ln[l-ze-'^^M 



(Al) 



ze 



-P^k 
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where the primed sum omits the k = term and the 
fugacity z is defined by 



z-e^^ 



< z < 1. 



The particle number is 






dz 



(A2) 



^^J_M<E1^_,^,^_V^A^ (A3) 



The critical point is where z — > 1, in which case 



v = vnQ + F{\,T) 



(A5) 



is an equation giving the condensate fraction no as a func- 
tion of T. The maximum value of T = Tg is when no = 0, 
and is the solution of the equation 



with A^o = 2:/(l — z). Dividing by N^ and replacing the 
sums over k by an integral over d^g gives 



V = vuq + F{z,T) , 



ofSCk 



+1 



d-'q 



(A4a) 



el^'^JT _ ^ 



u^Fil.To) 



(A6) 



for fixed value of i'. For i^ = 1, we have To/ J — 5.591 as 
expected [5]. 
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